setwd("/home/uec-00/yapingli/code/Allele_specific_methylation/")
for (e in commandArgs(TRUE)) {
  ta = strsplit(e,"=",fixed=TRUE)
  if(! is.na(ta[[1]][2])) {
		if(ta[[1]][1] == "k"){
			k<-ta[[1]][2]
			
		}
		if(ta[[1]][1] == "s"){
			s<-ta[[1]][2]
			}
	}
}

		fileName=paste("methylCGsRich_ASM_AllSnp_",k,"Merge_chr",s,"_table.txt",sep="")
		x<-read.table(fileName,header=F,sep="\t")
		x<-x[rowSums(x[,3:4])>=2 & rowSums(x[,5:6])>=2 & rowSums(x[,3:6])>=5 & rowSums(x[,3:6])<=50,]
		ASM<-array()
		i=1
		while(i<=length(x[,2])){
			temp<-x[i,1]
			y<-x[x[,1]==temp,]
			len<-length(y[,2])
			
			y<-cbind(y,rowSums(y[,3:6]))
			if(length(y[,3])>=15){
				y<-y[order(-y[,7]),]
				y<-y[1:15,]
			}
			minRange<-min(y[,2])
			maxRange<-max(y[,2])
			region=paste(minRange,maxRange,sep="~")
			A<-c(y[,3],y[,4])
			B<-c(y[,5],y[,6])
			
			p<-fisher.test(rbind(A,B),workspace = 2000000000)
			logP<-log2(p$p.value)
			methy_A=0
			methy_B=0
			for(j in c(1:length(y[,2]))){
				methy_A <- methy_A + y[j,3]/(y[j,3]+y[j,4])
				methy_B <- methy_B + y[j,5]/(y[j,5]+y[j,6])
			}
			methy_ratio=0
			ifelse(methy_B != 0,methy_ratio<-methy_A/methy_B,methy_ratio<-c(1000))

			thisRow<-cbind(temp,region,length(y[,2]),methy_ratio,p$p.value)
			ASM<-rbind(ASM,thisRow)

			i<-i+len
		}

		pAdjust<-p.adjust(as.numeric(ASM[,5]),method="BH")
		pAdjustLog<-log2(pAdjust)
		ASM<-cbind(ASM,pAdjust)
		outputName=paste("methylCGsRich_ASM_AllSnp_",k,"Merge_chr",s,"_pValue_ASMblock.txt",sep="")
		write.table(ASM, outputName,sep="\t",quote=F,col.names=F,row.names=F)


